
#ifndef __moab_project_soln_h__
#define __moab_project_soln_h__

#include "projection_errors.h"

// MOAB includes
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/ParallelComm.hpp"


#define PRINT_LAST_ERROR \
    if (MB_SUCCESS != result) {\
      std::string tmp_str;\
      std::cout << "Failure; message:" << std::endl;\
      moab->get_last_error(tmp_str);\
      std::cout << tmp_str << std::endl;\
      print_stacktrace(); \
      CHKERRABORT(PETSC_COMM_WORLD, result); }


// declare our functions

// int libmesh_project_solution(EquationSystems& src, EquationSystems& target, unsigned int dim=3, bool verbose=false, std::string output_file="");

/**
 * This is the exact solution functor.
 */
class MOABProjection : public ProjectionBase
{
public:
  MOABProjection(bool is_verbose) ;

  virtual ~MOABProjection() ;

  virtual int apply(std::string output_file="") ;

  virtual std::string type() ;

  void use_new_moab_method(bool) ;

private:

  int find_moab_libmesh_index_map(EquationSystems& es, moab::Core* moab, const moab::Range& verts, std::map<int,int>& moab_libmesh_map) ;

  moab::ErrorCode report_iface_ents(moab::Interface *mbImpl,
                              std::vector<moab::ParallelComm *> &pcs,
                              const bool print_results) ;

  /* Taken from mbcoupler_test.cpp */
  moab::ErrorCode test_interpolation(moab::Interface *mbImpl,
							               unsigned int dim,
                             std::string &interpTag,
                             std::string &gNormTag,
                             std::string &ssNormTag,
                             std::vector<const char *> &ssTagNames,
                             std::vector<const char *> &ssTagValues,
                             std::vector<moab::EntityHandle> &roots,
                             std::vector<moab::ParallelComm *> &pcs,
                             double &instant_time,
                             double &pointloc_time,
                             double &interp_time,
                             double &gnorm_time,
                             double &ssnorm_time,
                             double & toler,
                             bool verbose=false) ;

  int copy_from_to_moab_tag(moab::Core* moab, EquationSystems& es, moab::Tag tagh,
                            moab::Range& entities, bool to_moab,
                            std::map<int,int>& moab_libmesh_map) ;

  int copy_from_to_moab_tag(moab::Core* moab, EquationSystems& es, std::string tag_name,
                            moab::Range& entities, bool to_moab,
                            std::map<int,int>& moab_libmesh_map) ;

  class internal_node_compare
  {
    public:
      internal_node_compare(const double* base_vtx) : moab_vtx(base_vtx), tolerance(1e-7)
      { }

      bool operator() ( const Node* libmesh_node)
      {
        //return ((*libmesh_node-moab_node).size_sq() < tolerance*tolerance) ;
        return (fabs((*libmesh_node)(0)-moab_vtx[0]) < tolerance && fabs((*libmesh_node)(1)-moab_vtx[1]) < tolerance
              && fabs((*libmesh_node)(2)-moab_vtx[2]) < tolerance) ;
      }

    private:
      const double* moab_vtx;
      const double tolerance;
  };


  // static members
  static int mbconvert(int argc, const char* argv[]) ;

  static void reduceMax(double &v) ;

  // private data
  bool use_new_method;

};

inline
MOABProjection::MOABProjection(bool is_verbose) : ProjectionBase(is_verbose)
{
}

inline
MOABProjection::~MOABProjection()
{
}


inline
std::string MOABProjection::type()
{
  return "MOAB_Projection" ;
}


inline
void MOABProjection::use_new_moab_method(bool flg)
{
  use_new_method = flg;
}


#endif // __moab_project_soln_h__
